MODULE boots
  USE GLOBVAR
  USE RANDOM
  USE STATISTICS
  USE MATRIX
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE boot(momentboot)
    IMPLICIT NONE
    REAL(8), INTENT(INOUT):: momentboot(:,:) !weight(:,:) !initial guess
    REAL(8) :: a(N_ind,n_testc),momb(n_boot,n_mom),momb2(n_boot,n_mom2),wei(n_mom2,n_mom2)
    REAL(8) :: inter(n_ind,n_levelc,11,3),interc(n_ind,n_levelc,3)
    REAL(8),ALLOCATABLE:: AVSE(:),cor(:,:)
    INTEGER :: duc(n_ind,n_levelc),du(n_levelc),ind(N_ind,n_levelc)
    INTEGER :: i,j,k,t,m,n,l,indexc(n_testc,2), artc(n_ind,n_levelc,11),mm(n_ind),ii,jj,tt,kk,mmm 
   
    
    !!!!! start to calculate data moments!!!!!
    IF (myid==0) THEN
    !!! language tasks!!!!
   DO ii=1,n_boot     
       jj=ii
        CALL SET_SEED(1438914432,234621+jj,1377777,123727)
       DO i=1,N_ind
           mm(i)=INT(SAMPLE_UNIFORM(1.0D0,DBLE(N_IND+1)))
           mmboot(ii,i)=mm(i)
       END DO
       
       
        !!! cognitive tasks!!!!
    indexc=0
    DO i=1,N_ind
      DO j=1,N_testc
         IF (c_test(mm(i),j)>-90) THEN
             indexc(j,1)=indexc(j,1)+1
             IF (c_test(mm(i),j)>0.5) indexc(j,2)=indexc(j,2)+1
         END IF
      END DO
    END DO
    
    DO i=1,N_testc
        IF (indexc(i,1)>0.5) THEN
            momb(ii,i)=DBLE(indexc(i,2))/DBLE(indexc(i,1)) ! PASSING RATE FOR EACH Cognitive TASK
        ELSE 
            momb(ii,i)=-99
        END IF
    END DO
    
    !!!! generate artf variables !!!!
    artc=-1
    DO i=1,N_ind
        DO j=1,N_testc
            IF (c_test(mm(i),j)>-90) THEN
                k=levelc(j)
                m=repc(mm(i),j)
                artc(i,k,m)=c_test(mm(i),j)
                inter(i,k,m,1:3)=x_testc(mm(i),j,1:3)
            END IF
        END DO
    END DO
    
    !!!! passing rate at each level!!!!

    DO j=N_testc+1,N_testc+N_levelc
       m=0
       n=0 
       DO i=1,N_ind
           DO k=1,11
               IF (artc(i,j-N_testc,k)>-0.5) THEN
                   m=m+1
                   IF (artc(i,j-N_testc,k)==1) n=n+1
               END IF
           END DO 
       END DO
       IF (m>0) THEN
           momb(ii,j)=DBLE(n)/DBLE(m)
       ELSE
           momb(ii,j)=-99
       END IF
    END DO
    
    !!!! correlation across levels !!!!
    
    mmm=N_testc+N_levelc
    
     l=0
    DO j=1,N_levelc-2
        IF (j<=2 .or. (j>4 .and. j<=8)) THEN
            tt=5
        ELSE IF (j>2 .and. j<5) THEN
            tt=2
        ELSE IF (j==9) THEN
            tt=2
        ELSE IF (j==10) THEN
            tt=4
        ELSE IF (j==11) THEN
            tt=5
        END IF
        
        IF (j<2 .or. (j>3 .and. j<8)) THEN
            kk=5
        ELSE IF (j>1 .and. j<4) THEN
            kk=2
        ELSE IF (j==8) THEN
            kk=2
        ELSE IF (j==9) THEN
            kk=4
        ELSE IF (j==10) THEN
            kk=5
        ELSE IF (j==11) THEN
            kk=3
        END IF
        
         DO t=1,tt
           DO k=1,kk
             m=0
             n=0 
             l=l+1
             DO i=1,N_ind 
                IF (artc(i,j,t)==1 .and. artc(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                    m=m+1
                    IF (artc(i,j,t)==1 .and. artc(i,j+1,k)==1) n=n+1
                END IF
             END DO
             
             IF (m>0) THEN
                momb(ii,mmm+L)=DBLE(n)/DBLE(m)
             ELSE
                momb(ii,mmm+L)=-99
             END IF
           END DO
         END DO
    END DO
    
    mmm=N_testc+N_levelc+177
    
    l=0
    DO j=1,N_levelc-1
        IF (j<=2 .or. (j>4 .and. j<=8)) THEN
            tt=4
            kk=5
        ELSE IF (j>2 .and. j<5) THEN
            tt=1
            kk=2
        ELSE IF (j==9) THEN
            tt=1
            kk=2
        ELSE IF (j==10) THEN
            tt=3
            kk=4
        ELSE IF (j==11) THEN
            tt=4
            kk=5
        ELSE IF (j==12) THEN
            tt=2
            kk=3
        END IF
        
         DO t=1,tt
           DO k=t+1,kk
             m=0
             n=0 
             l=l+1
             DO i=1,N_ind 
                IF (artc(i,j,t)==1 .and. artc(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                    m=m+1
                    IF (artc(i,j,t)==1 .and. artc(i,j,k)==1) n=n+1
                END IF
             END DO
             
             IF (m>0) THEN
                momb(ii,mmm+L)=DBLE(n)/DBLE(m)
             ELSE
                momb(ii,mmm+L)=-99
             END IF
           END DO
         END DO
    END DO
    
     mmm=N_testc+N_levelc+177+82
    
    l=0
    DO j=1,N_levelc-3
        IF (j==1) THEN
            tt=5
            kk=2
        ELSE IF (j==2) THEN
            tt=5
            kk=2
        ELSE IF (j==3) THEN
            tt=2
            kk=5
        ELSE IF (j==4) THEN
            tt=2
            kk=5
        ELSE IF (j==5) THEN
            tt=5
            kk=5
        ELSE IF (j==6) THEN
            tt=5
            kk=5
        ELSE IF (j==7) THEN
            tt=5
            kk=2
        ELSE IF (j==8) THEN
            tt=5
            kk=4
        ELSE IF (j==9) THEN
            tt=2
            kk=5
        ELSE IF (j==10) THEN
            tt=4
            kk=3
        END IF
        
        
      DO t=1,tt
        DO k=1,kk
          m=0
          n=0 
          l=l+1
          DO i=1,N_ind 
             IF (artc(i,j,t)==1 .and. artc(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artc(i,j,t)==1 .and. artc(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             momb(ii,mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             momb(ii,mmm+l)=-99
          END IF
        END DO
      END DO
    END DO  
    
     mmm=N_testc+N_levelc+177+82+142
    
    
    indexc=0
    DO i=1,N_ind
      DO j=1,N_testc
          a(i,j)=firstc(mm(i),2)-monthc(j)-DBLE(timec(j))/4.0d0
         IF (c_test(mm(i),j)>-90 .and. a(i,j)>-1.25d0) THEN
             indexc(j,1)=indexc(j,1)+1
             IF (c_test(mm(i),j)>0.5) indexc(j,2)=indexc(j,2)+1
         END IF
      END DO
    END DO
    
    DO i=1,N_testc
        IF (indexc(i,1)>0.5) THEN
            momb(ii,mmm+i)=DBLE(indexc(i,2))/DBLE(indexc(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            momb(ii,mmm+i)=-99
        END IF
    END DO
    
    mmm=2*N_testc+N_levelc+177+82+142
    
    indexc=0
    DO i=1,N_ind
      DO j=1,N_testc
          a(i,j)=DBLE(firstc(mm(i),2)-monthc(j)-DBLE(timec(j))/4.0d0)

         IF (c_test(mm(i),j)>-90 .and. a(i,j)<-1.25d0) THEN
             indexc(j,1)=indexc(j,1)+1
             IF (c_test(mm(i),j)>0.5) indexc(j,2)=indexc(j,2)+1
         END IF
         
      END DO
    END DO

    
    
    DO i=1,N_testc
        IF (indexc(i,1)>0.5) THEN
            momb(ii,mmm+i)=DBLE(indexc(i,2))/DBLE(indexc(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            momb(ii,mmm+i)=-99
        END IF
    END DO
    
    mmm=3*N_testc+N_levelc+177+82+142
    
    l=0
    DO j=1,N_levelc-1
      DO t=1,5
          m=0
          n=0
          l=l+1
          DO i=1,N_ind 
             IF (artc(i,j,t)>-0.5) THEN   
                 m=m+1
                 IF (artc(i,j,t)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             momb(ii,mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             momb(ii,mmm+l)=-99
          END IF
      END DO
    END DO
    mmm=3*N_testc+N_levelc+177+82+142+5*(N_levelc-1)  
    
    DUC=-99
    ind=-99
    DO i=1,N_ind
        DO k=1,N_levelc
            DO j=1,11
                IF (artc(i,k,j)==1 .and. j==1) THEN
                    duc(i,k)=1
                    ind(i,k)=1
                ELSE IF (artc(i,k,j)==0 .and. j==1) THEN
                    duc(i,k)=j
                    ind(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF (artc(i,k,j-1)==0 .and. artc(i,k,j)==0 .and. ind(i,k)==0) THEN
                       duc(i,k)=j
                       ind(i,k)=0
                       IF (j<11) THEN
                           IF (artc(i,k,j+1)==-1) THEN
                           duc(i,k)=j+1
                           ind(i,k)=1
                           END IF
                       END IF                      
                    ELSE IF (artc(i,k,j-1)==0 .and. artc(i,k,j)==0 .and. ind(i,k)==1) THEN
                        duc(i,k)=duc(i,k)
                        ind(i,k)=1
                    ELSE IF (artc(i,k,j-1)==0 .and. artc(i,k,j)==1 .and. ind(i,k)==0) THEN
                        duc(i,k)=j
                        ind(i,k)=1
                    ELSE IF (artc(i,k,j-1)==0 .and. artc(i,k,j)==1 .and. ind(i,k)==1) THEN
                        duc(i,k)=duc(i,k)
                        ind(i,k)=1
                    END IF
                END IF
            END DO           
        END DO
    END DO
    
    DO k=1,n_levelc
        m=0
        n=0
        DO i=1,n_ind
            IF (duc(i,k)>-90) THEN
               m=m+1
               n=n+duc(i,k)
            END IF
        END DO
        IF (m>0) THEN
          momb(ii,mmm+k)=DBLE(n)/DBLE(m)
        ELSE 
           momb(ii,mmm+k)=-99
        END IF
    END DO
    
    mmm=3*N_testc+N_levelc+177+82+142+5*(N_levelc-1)+n_levelc
    
    du=0
    DO k=1,n_levelc
        DO i=1,N_ind
            IF (duc(i,k)>-90) THEN
                du(k)=du(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levelc
        allocate(avse(du(k)))
        m=0
        DO i=1,n_ind
            IF (duc(i,k)>-90) THEN
                m=m+1
                avse(m)=duc(i,k)
            END IF
        END DO
        momb(ii,mmm+k)=STDEV_V(AVSE)
        DEALLOCATE(avse)
    END DO
    mmm=3*N_testc+N_levelc+177+82+142+5*(N_levelc-1)+2*n_levelc
    

    IND=0
    interc=0.0d0
    DO k=1,N_levelc
        DO i=1,N_ind
          DO j=1,11
              IF (inter(i,k,j,1)>-90) THEN
                  ind(i,k)=ind(i,k)+1
                  interc(i,k,1)=interc(i,k,1)+inter(i,k,j,1)
                  interc(i,k,2)=interc(i,k,2)+inter(i,k,j,2)
                  interc(i,k,3)=interc(i,k,3)+inter(i,k,j,3)
              END IF
          END DO
          IF (ind(i,k)>0) THEN
               interc(i,k,1)=DBLE(interc(i,k,1))/DBLE(ind(i,k))
               interc(i,k,2)=DBLE(interc(i,k,2))/DBLE(ind(i,k))
               interc(i,k,3)=DBLE(interc(i,k,3))/DBLE(ind(i,k))
          END IF
        END DO
    END DO
    
    
    du=0
    DO k=1,n_levelc
        DO i=1,N_ind
            IF (ind(i,k)>0 .AND. duc(i,k)>-90) THEN
            du(k)=du(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levelc
        allocate(cor(du(k),4))
        m=0
        DO i=1,n_ind
             IF (ind(i,k)>0 .AND. duc(i,k)>-90) THEN
                m=m+1
                cor(m,1)=interc(i,k,1)
                cor(m,2)=interc(i,k,2)
                cor(m,3)=interc(i,k,3)
                cor(m,4)=duc(i,k)
             END IF
        END DO
        IF (k<n_levelc) THEN
        momb(ii,mmm+3*(k-1)+1)=Correlation(cor(:,1),cor(:,4))
        momb(ii,mmm+3*(k-1)+2)=Correlation(cor(:,2),cor(:,4))
        momb(ii,mmm+3*(k-1)+3)=Correlation(cor(:,3),cor(:,4))
        ELSE 
            momb(ii,mmm+3*(k-1)+1)=-99
            momb(ii,mmm+3*(k-1)+2)=-99
            momb(ii,mmm+3*(k-1)+3)=-99
        END IF
        DEALLOCATE(cor)
    END DO
    mmm=3*N_testc+N_levelc+177+82+142+5*(N_levelc-1)+5*n_levelc
    
      
    
   END DO !ii
   
   momentboot=momb
   
       
           !
           !wei=Variance_Covariance(momb2)
           !!weight=Matrix_Inverse_symmetric(wei)
           !weight=Matrix_Inverse(wei)
           !OPEN(12,file="wei.out")
           !  DO j=1,N_mom2
           !   WRITE(12,'(2000F32.12)') (wei(j,i),i=1,n_mom2)
           !  END DO
           !CLOSE(12) 
           !OPEN(12,file="weight.out")
           !  DO j=1,N_mom2
           !   WRITE(12,'(2000F32.12)') (weight(j,i),i=1,n_mom2)
           !  END DO
           !CLOSE(12)  
           !OPEN(12,file="momm.out")
           !  DO i=1,N_mom
           !   WRITE(12,'(2000F32.12)') (momb(j,i),j=1,n_boot)
           !  END DO
           !CLOSE(12) 
              OPEN(12,file="mmboot.out")
             DO i=1,N_ind
              WRITE(12,'(2000I4)') (mmboot(j,i),j=1,n_boot)
             END DO
           CLOSE(12) 
      
    
    
    END IF  
        
    
    
    
    
  END SUBROUTINE boot
  
END MODULE boots